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Abstract 

We consider instability of a two layered solid body of a denser mate- 
rial on top of a lighter one. This problem is widely known to geoscientist 
in sediment-salt migration as salt diapirism. In the literature, this prob- 
lem has often been treated as Raleigh- Taylor instability in viscous fluids 
instead of solid bodies. In this paper, we propose a successive linear ap- 
proximation method for large deformation in viscoelastic solids as a model 
for salt migration. 

Keywords: Large deformation, viscoelastic solids, successive linear ap- 
proximation, boundary value problem, incremental method, numerical 
simulation. 

1 Introduction 

The behavior for large deformation is characterized by some nonlinear constitu- 
tive relations, which leads to a system of nonlinear partial differential equations. 
To solve boundary value problems of this nonlinear system for large deforma- 
tion, we propose a method in a successively updated referential formulation of 
linear approximation. 

The method is based on the well-known problem of small deformation su- 
perposed on finite deformation in the literature [1]. Roughly speaking, at each 
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step, the constitutive functions are calculated at the present state of deforma- 
tion which will be regarded as the reference configuration for the next state, and 
assuming the deformation to the next state is small, the constitutive functions 
can be linearized. In this manner, it becomes a linear problem from one state 
to the next state successively with small deformations. 

Numerical simulations for elastic bodies with this method have been success- 
fully compared with some exact solutions for large deformations (2j [3] . In this 
paper, we consider the instability of a two-layered solid body of a denser elastic 
material on top of a lighter viscoelastic one. This problem is widely known 
to geophysicists in sediment-salt migration as salt diapirism. Due to techni- 
cal difficulties in dealing with large deformations in solid bodies, the problem 
has often been treated as Raleigh- Taylor instability in viscous fluids (see for 
example, [HE]) instead of solid bodies. 

We will show in our numerical simulation with finite element method the 
formation of diapirs with viscoelastic solid models. The mesh of the calculation 
domain is updated as the body deforms at every step, however, up to very large 
deformation we have found that re-meshing is almost unnecessary. 

Remark 1: 

For numerical computation of large deformations in finite elasticity, to obviate 
the difficulty of handling nonlinearities, incremental methods are widely dis- 
cussed, for example, in the books by Ciarlet [BJ, Oden [7J, and Ogden jjj. These 
methods consist in letting the boundary forces vary by small increments from 
zero to the given ones and to compute corresponding approximate solutions by 
successive linearization. The general idea seems to be of purely mathematical 
concern regarding linearization between successive boundary value problems. 
The problems are usually formulated with domains in the initial reference con- 
figuration, i.e., in (total) Lagrangian formulation. This is the essential difference 
from the method proposed in this paper with "updated" Lagrangian formula- 
tion, namely, the reference configuration is updated to the present state at each 
step. However, we refrain from calling the proposed method as an updated La- 
grangian formulation (UL) , because UL formulation is widely known in numer- 
ical methods albeit frequently with quite different tenets (see for example [9]) 
in which the basic equations and numerical treatment of boundary conditions 
are quite different. 

The present approach has an additional advantage that the assignment of 
boundary data is simpler and straightforward. They are prescribed on the 
present state and the difference between prescribing the corresponding bound- 
ary data on the surface of the consecutive states is of higher order which is 
insignificant in the linear approximation. Therefore, unlike the usual incremen- 
tal method in Lagrangian formulation, which needs to prescribe the incremental 
boundary data on the surface of the initial reference state with values depend- 
ing on the deformed geometry of the body, we are totally free from any such 
complicated concerns (see [HOI])- 
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Remark 2: 

For large deformations, systems of governing partial different equations are gen- 
erally nonlinear. The boundary value problems are usually formulated in ref- 
erential (Lagrangian) or spatial (Eulcrian) coordinates, and numerical solution 
of nonlinear systems of algebraic equations by Newton's type methods. In this 
paper, we formulate boundary value problems in the coordinates relative to the 
configuration at the present time. In other words, we shall describe deformations 
relative to the current configuration, instead of a fixed reference configuration 
in the Lagrangian formulation. Note that this is not an Eulerian formulation 
either. The Eulerian coordinates are fixed coordinates in space. We shall call 
this a relative-descriptional formulation [3] . 

The advantage of using the relative-descriptional formulation is that we can 
consider a small time step from the present state, so that the constitutive func- 
tions can be linearized relative to the present state, and the equation of motion 
becomes linear in relative displacement. When the present state proceeds in 
time, a nonlinear finite deformation can be treated as a sequence of small defor- 
mations in the same manner as the usual Euler's method for solving differential 
equations, i.e., successively at each state, the tangent is calculated and used to 
extrapolate the neighboring state. 

The essential idea of the proposed method is based on the approach of small 
deformation superposed on large deformation [I]. Although the small- on-large 
idea is very well-known, to my knowledge, numerical implementations of Euler's 
type successive approximation based on this idea have not been explored. 

On the other hand, computational literature based on the small-on-large 
idea, typically the methods introduced in [10] , [IT] , and [12] are very well docu- 
mented. Such methods, in which the nonlinear algebraic systems resulting from 
the variational formulation of the nonlinear boundary value problem are solved 
by Newton's method, employ the tangent operator obtained from small-on-large 
linearization of the constitutive functions and boundary conditions. Such meth- 
ods are quite different from the one proposed in this paper concerning the use 
of the small-on-large approach. 

2 The present configuration 

Let Kg be a reference configuration of the body B, and Kt be its deformed 
configuration at the present time t, Let 



x = X (X,t) 



X e k {B), 



and 




£ = x(X,t) :=£(*, t) Gk t (B), 



x = x(X,t) e K t (B). 



3 



It can also be regarded as the relative deformation at time r with respect to the 
present configuration at time t denoted as £(x,t). We also define the relative 
displacement vector as 

u(x, t) = £(x, t) — x. 

Note that 

Vx«x,r) = V x ( X (X,t)) V^x^t))- 1 = F(X,r)F(X,t)-\ 
hence, we have 

V x w(x,r) = F(X,r)F(X,t)- 1 - /, 

or simply as 

F(r) = (I + H(r))F(t), (I) 
where I is the identity tensor and 

H(x,t) = V x u{x,t) (2) 

is the relative displacement gradient at time r with respect to the present con- 
figuration Kt (emphasize, not kq). Moreover, by taking the time derivative with 
respect to r, it gives 

F(t) = H(r)F(t). (3) 

In these expressions and hereafter, we shall often denote a function F as F(t) 
to emphasize its value at time t when its spatial variable is self-evident. 

In summary, we can represent the deformation and deformation gradient 
schematically in the following diagram: 




x e K t (B) 

- X + u 



3 Linearized constitutive equations 

Let k be the preferred reference configuration of a viscoelastic body B, and let 
the Cauchy stress T(X, t) be given by the constitutive equation in the configu- 
ration Kq, 

T(X,t) = T(F(X,t),F(X,t)). (4) 

For large deformations, the constitutive function T is generally a nonlinear 
function of the deformation gradient F. 

We shall regard the present configuration Kt as an updated reference config- 
uration, and consider a small deformation relative to the present state Kt(B) at 
time t = t + At. In other words, we shall assume that the relative displacement 
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gradient H is small, \H\ <C 1, so that we can linearize the constitutive equation 
(4) at time r relative to the updated reference configuration at time t, namely, 

T(r) = T(F(r),F(r)) =T(F(t),0) 

+ d F T(F(t),0)[F(r) - F(t)] + d p T(F(t),0)[F(r)] + o(2), 

or by the use of (3), 

T(T)^T e (t)+d F r(F(t),0)[H(T)F(t)]+d p T(F(t),0)[H(T)F(t)]+o(2). 

where 

T e (t)=T(F(t),0) 

is the elastic Cauchy stress at the present time t and o(2) represents higher 
order terms in the small displacement gradient \H\. 

The linearized constitutive equation can now be written as 

T(r) = T e (t) + L(F(t))[H(r)} + M(F(t))[H(r)], (5) 

where 

L(F)[H]:=d F T(F,0)[HF], M(F)[H] := d F T(F,0)[HF], (6) 

define the fourth order elasticity tensor L(F) and viscosity tensor M(F) relative 
to the present configuration Kt- 

The above general definition of the elasticity and the viscosity tensors for any 
constitutive class of viscoelastic materials T — T(F, F), relative to the updated 
present configuration, will be explicitly determined in the following sections for 
a particular class, namely a Mooney-Rivlin type materials. 



3.1 Compressible and nearly incompressible bodies 

For a viscoelastic body, without loss of generality, the constitutive equation (4) 
relative to the preferred reference configuration n can be written as 

T = T(F,F) = -p(F)I + f(F,F). (7) 

However, for an incompressible body, the pressure p depends also on the bound- 
ary conditions, and because it cannot be determined from the deformation of 
the body alone, it is called an indeterminate pressure, which is an indepen- 
dent variable in addition to the displacement vector variable for boundary value 
problems. 

For compressible bodies, we shall assume that the pressure depend on the 
deformation gradient only through the determinant, or by the use of the mass 
balance, depend only on the mass density, 

p = p(dctF) =p(p), p= f° , 
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where po is the mass density in the reference configuration k - We have 

p(r) =p (dctF(r)- 1 -detF(i)- 1 ) =p(i)(det(F(r)- 1 F(i)) -1) 

= p(i)(det(J + H(t))- 1 - 1) = -p(t) trJT(r) + o(2), 

in which the relation (1) has been used. 
Therefore, it follows that 

P(r)-P(t) = (^) t (p(T)-p(t))+o(2) = -(p|) t trH(r) +o(2), 

or 

p(T)=p(t)-ptvH(r) + o(2), (8) 
where /3 := (p is a material parameter evaluated at the present time t. 
From (7) and (1), let 

C(F(t))[H(r)} :=d F f(F(t),0)[F(r)-F(t)]=d F f(F(t),0)[H(r)F(t)}, 

then from (6)i, the elasticity tensor becomes 

L(F)[H]=f3(tvH)I + C(F)[H}. (9) 

We call a body nearly incompressible if its density is nearly insensitive to the 
change of pressure. Therefore, if we regard the density as a function of pressure, 
p = p(p), then its derivative with respect to the pressure is nearly zero. In other 
words, for nearly incompressible bodies, we shall assume that (5 is a material 
parameter much greater than 1, 

/3»1. (10) 

Note that for compressible or nearly incompressible body, the elasticity ten- 
sor does not contain the pressure explicitly. It is only a function of the defor- 
mation gradient and the material parameter (i at the present time t. 

3.2 A viscoelastic material model 

As an example, we shall consider a nearly incompressible viscoelastic material 
with the following constitutive equation, 

T = -pI + f(F,F), 

f(F,F) = Sl B + s 2 B- 1 

+ A(tr D)I + + p 2 (DB + BD) + p^DB- 1 + B^D), 

where B = FF T is the left Cauchy-Green strain tensor and D is the rate 
of strain tensor. The material parameters s± through /13 are assumed to be 
constants. This material model will be referred to as a Mooney-Rivlin type 
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isotropic viscoelastic solid, since, by assuming the relevant material parameters 
to be constant, the elastic part of the constitutive equation represents the well- 
known Mooney-Rivlin isotropic elastic solid in nonlinear elasticity (see [TSim] ). 
Note that this constitutive equation contains only linear terms in D (see also 
|15jV Therefore, it would be appropriate for large deformations with small strain 
rate. 

After taking the gradients of T(F, F) with respect to F at (F, 0), we have 
C(F)[H] = Sl (HB + BH T ) - s 2 (B- 1 H + H T B^ 1 ), 
which by (9), gives the elasticity tensor 

L(F)[H] = -P(ttH)I + sx(HB + BH T ) - s 2 (B~ 1 H + H T '.B^ 1 ), (12) 
and with respect to F, we obtain 

M(F)[H] = A(tr H)I + M (H + H T ) + (H + H T )M , (13) 

where 

M o ■= 2 {(J-iI + ^B + ^B- 1 ). 

From (5), we have 

T(r) = T e {t) + L(F(t))[H(r)} + M (F '(t))[H (r)] . (14) 

Furthermore, the (first) Piola-Kirchhoff stress tensor at time r relative to the 
present configuration at time t, denoted by T*(r), is given by 

%(t) = det(7 + H) T(t)(I + H)~ T 

= det(J + H){T e {t) + L(F)[H] + M(F)[H])(I + H)- T 

= (I + tr H) (T e (t) + L(F) [H] + M(F) [H]) (I — H T ) + o(2) 

= T e {t) + (tr H)T e (t) - T e (t)H T + L(F){H] + M(F)[H] + o(2). 

We can write the linearized Piola-Kirchhoff stress as 

T t (r) = T e (t) + (tr H(r))T e (t) - T e {t)H{r) T 
+ L(F(t))[H(r)}+M(F(t))[H(r)}. 

Note that when r — ► t, H — > 0, and hence T t (r) — > T(t), therefore, the Piola- 
Kirchhoff stress, becomes the Cauchy stress at the present time t. 
We can also write 

T t (T)=T e (t)+K(F(t),T e (t))[H(T)]+M(F(t))[H( T )], (15) 

where the Piola-Kirchhoff elasticity tensor is defined as 

K(F, T e )[H] := (tr H)T e - T e H T + L(F)[H\. (16) 

In numerical examples presented later, the material is assumed to be of 
Mooney-Rivlin type and these relations will be used. 
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4 Boundary value problem 



Let fl = Kt(B) C M 3 be the region occupied by the body at the present time t, 
and d£l = Ti U r 2 , and n K be the exterior unit normal to dVt. We can now 
consider the updated referential formulation of the boundary value problem. 

At time r > t, we shall consider the boundary value problem in the La- 
grangian formulation, with the present state at time t as the updated reference 
configuration, given by 

-divT t (r)=p(t) ff (T), infi, 
T t (r)n K = f(r), on T 1; 
u(r) • n K = 0, on T 2 , 
, T t (r)n K x n K = 0, on T 2 , 

where Tt(x, r) is the Piola-Kirchhoff stress at time t with respect to the state 
at the present time t. The divergence operator is relative to the coordinate 
system (x) in the present configuration. The body is subjected to the body 
force g(x,r) and the surface traction f(x,r). 

From the definition (2), we also have the initial condition for the displace- 
ment vector u(x,t): 

u(x,t)=0, VxeK t (B). (18) 

Since Tt{r)n R x n K = implies that the surface traction T t (r)n K is in the 
direction of the normal, the last boundary condition in (17) states that the 
tangential component of the surface traction T t (r)n K vanishes on r 2 . In other 
words, the boundary T 2 is a roller-supported boundary. 

4.1 Linearized boundary value problem 

We shall assume that at the present time t, the deformation gradient F with 
respect to the preferred reference configuration kq and the Cauchy stress T 
are known, and that r = t + At with small enough At, then from (15), the 
equilibrium equation in (17) can be written as 

-div(K(F(t),T e (t))[Vu(T)} + M(F(t))[Vu(T)}) = divT e (t) + p(t)g{r), (19) 

which is a linear partial differential equation for the displacement vector u(x, r). 
The right hand side is a known function. 

The idea of formulating the boundary value problem in the form (17) for 
small time increment is similar to the theory of small deformations superposed 
on finite deformations (see [U [14] ) . 

4.2 Variational formulation 

The boundary value problem (17) can be formulated as a variational problem. 
Let us consider the Sobolev space of vector valued functions on f2, 

H 1 ^) = {v : fl -> M 3 | v, Vv e L 2 (n)} 
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and the subspace 

V = {v e H\n) | v ■ n K = on T 2 }. 

Taking the inner product of the equation (19) with a vector w 6 V and 
integrating over the domain fi, we obtain, after integration by parts, 

K[Vu(t)] + M[V«(r)]) • Vwdv 

= / p{t)g{r) ■ w dv - / T e (t) -Vwdv + I T t {r)n K ■ w da, 
Jn Jn Jon 

where the relation (15) is used and dot ( • ) represents the inner product between 
two vectors as well as the inner product between two second order tensors, i.e., 
A - B = tv(AB T ). 

Since w G V, w ■ n K = 0, by the boundary conditions the surface integral 
on the right hand side becomes 

/ T t (r)n K ■ w da = / f(j)-wda. 
Therefore, if we define the following bilinear forms 

C(w,u(t)) := / K[Vu(t)} -Vwdv, 

Jn (20) 

M(w,u(t)) := / M[Vu(t] ■ Vwdv, 
Jn 

and the linear form, 

V{w):= I p{t)g{r) ■ w dv + / f(r)-wda- / T e {t)-Vwdv, (21) 
Jn Jr ± Jn 

then the variational problem is to find the solution vector u(t) & V such that 
C(w, u(t)) +M(w, u(t)) = V{w) y W eV. (22) 
Note that from the initial condition (18), u(x,t) = 0, we can approximate 

and hence restate the variational problem as: Find the solution vector u(t) e V 
such that 

K(w,u(t)) =V(w) VweV, (23) 

where 

K.(w, u) := C(w, u) + ~^.M( W ' u )- 

The variational problem depends on the elasticity and viscosity tensor at the 
updated present state. Mathematical analysis of requirements at the present 
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state for existence and uniqueness of solution will be presented in a forthcom- 
ing paper. In general, non-existence or non-uniqueness may occurred if such 
requirements are not fulfilled. 

For numerical solutions of the variational equation, finite element method 
will be used. 

5 Successive linear approximation 

Recall the Euler method of solving differential equation, say y — f(t), that for 
a discrete time axis, • • • < t n -\ < t n < t n +i <■•■■> an d y(t n ) = Vn, the solution 
curve can be constructed by y„+i = y n + f(t n )At, where f(t n ) is the tangent 
of the solution curve at t n . We can use a similar strategy for solving problems 
of large deformation, by solving linear variational problem stated in (23). 

We consider a discrete time axis, to < ■ ■ ■ < t n < t n+ \ < • • • with small 
enough time increment At. Let Kt n be the configuration of the body at the 
instant t n and 

x n = x{X,t n ) G K tn (B) for X G K r (B), n r = n to . 

Let the elastic Cauchy stress T e (x n , t n ) and the deformation gradient F(x n ,t n ) 
relative to the preferred configuration k at the present time t = t n be known. 
The boundary value problem (17) at the instant r = t n+ \ with the body force 
b(x n , t n+ i) and the surface traction f(x ni t n+1 ) in updated referential formula- 
tion with respect to the present configuration n tn , can now be solved numerically 
as a problem in linear elasticity for the relative displacement field u(x ni t n+1 ) 
from the present state at t n . 

After solving the problem (17) at t n , the configuration K tn+1 can be regarded 
as the reference configuration at the updated present time t n+ i from the dis- 
placement field, i.e., 

Xn+l = X(X,t n +l) = X n + u(x n ,t n+ i), 

while the deformation gradient 

F(x n+1 ,t n+1 ) = (I + Vu(x n ,t n+1 ))F(x n ,t n ) 

and the elastic Cauchy stress T e (x n+ i, t n+ i) can be calculated from (14) at 
t n+ i so that the updated referential formulation of the problem in the form 
(17), with the body force b{x n+ i, t n+ 2) and the surface traction f{x n+ \, t n+ 2), 
can proceed again from the updated referential configuration at t n+ \. This 
numerical procedure will be referred to as the successively updated referential 
formulation of linear approximation, or simply as the method of Successive 
Linear Approximation (SLA). The updating process can be represented in the 
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following schematic diagram: 



BVP _ 

/n\ _ fintcly K t — Kt n k t ~ K t n+1 

r[ ' deforms K tn (B) " + « tn+1 (B) 

update JJ- 



BVP 

«tn+l( fi ) K t n+2 (B) 

update JJ- 

«t = K t n + 2 



where BVP stands for the linear boundary value problem formulated in (17) to 
be solved successively after each update as shown in the above scheme, and the 
process starts with n t — K to = n r which is the initial reference configuration. 

5.1 Incremental loadings 

Note that the force term V(w) in the variational equation (23) defined in (21) 
is, in fact, a small quantity of the order of the time increment At. Therefore, 
the method of SLA can be regarded as an incremental method. Here we shall 
emphasize the incremental features of the boundary value problem (17). 

For convenience, we shall denote the time dependence on t n as subindex n 
or simply as n. For example, we write the Cauchy stress T(t n ) = T n and the 
elastic Cauchy stress T e (t n ) — T e (n), and we have 

T n = T{F n ,F n ) = T(F n , 0) + 8pT(F ni 0)[F n ] = T e (n) + dpT(F n ,0)[F n ] . (24) 

Now, we can rewrite the variational equation (23) as, 

JC(w,u n+1 ) =V(w,n) VweV, (25) 

where, from (21), the force term is given by 

V{w,n)= I p n g n+1 -wdv+ / f n+1 -wda- / T e (n) -Vwdv. (26) 
Jn JTr Jn 

Note that by (24), it follows that 

/ T e (n) -Vwdv= / T n -Vwdv- I d p T(F n , 0)[F n ] ■ Vw dv, (27) 
Jn Jn Jn 



and 



/ T n ■ Vw dv = T n n K ■ w dv — / div T n -wdv. 
Jn Jon Jn 
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On the other hand, noting that T tri (t n ) — T n is the Cauchy stress at t n , from 
(17) we have the following boundary value problem at the present time t n , 

-divT„ = p n g n , in Q, 
T n n K = f n , on Ti, 
^ T n n K x n K = 0, on T 2 , 

which leads to 

/ T n ■ Vw dv = f n wdv+ p n g n ■ w dv. (28) 
Jn Jr ± Jn 

Combing (26), (27), and (28), we obtain 

V(w,n) =h(w,n) + I 2 (w,n)+ h{w,n), (29) 
which contains three types of loading for the variational equation (25), namely, 

h(w,n) = / p n (g n+1 -g n )-wdv, 
Jn 

h(w,n) = / (f n+1 - /„) -wda, (30) 

I 3 (w,n)= [ dpT(F n ,0)[F n ]-Vwdv. 
Jn 

The integral I± represents the incremental body force and Iq. is the incremental 
surface traction between time steps t n and t n+ \. The third one -Z3 is due to the 
viscous effect of the material body. 



5.2 Incremental approximation for large deformations 

If we assume that the functions g(x,t) and f(x,t) are smooth in t, then their 
increments, <7 n+1 — g n and /„ +1 — /„, are of the order of At = t n+ \ — t n , and 
since the variational equation is linear, the solution vector u(x, t n +i) is also of 
the same order. For elastic material bodies, these are the two possible types of 
incremental loading. 

Starting from an initial solution and applying the method of SLA with proper 
loading conditions at each time step, the problem of large deformation can 
be obtained. Numerical examples employing the SLA method for incremental 
surface traction of pure shear and of gradually bending a rectangular block 
into a circular section have been considered in [2l [3] for Mooney Rivlin elastic 
materials. 

For viscoelastic material bodies, there is a possible loading due to the inte- 
gral I3. To see this, we shall consider the case without surface traction / = 
and time- independent body force g so that I\ = 1% =0, hence the integral ^3 
is the only possible incremental loading. 
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To begin with, for n = 0, assume that at the initial time to, we have a static 
equilibrium solution so that we have the initial conditions: Fq = I and Fq = 0. 
Therefore, from (25), (29) and (30)3, we have 

K(w,ui) = h(w,Q) = 0, VweV, 

which implies that the solution vector Ui(x) = 0. Consequently, iii — -htiii = 
and Hi = Vtti = 0. Then, from (3), F\ = HiFq = 0, which, in turn, leads to 

fc(t»,ua) = J3(u7,l) = 0, VweV, 

and implies that «2 must also vanish and so on. In other words, if the initial 
solution is an equilibrium solution, the solution remains valid for all time. How- 
ever, this conclusion may not be true since the initial solution may not be a 
stable equilibrium solution in general. 

Therefore, in order to study the stability, a small perturbation of the initial 
solution is needed so that Ui ^ 0, and hence Is(w,l) ^ 0, to trigger the 
successive evolution of deformations. Two such examples are considered in the 
following sections. 



6 Salt migration 



As an example of large deformation, we shall consider a body consisting of two 
different layers initially, with the mass density of the upper layer, the overburden 
sediment, greater than that of the bottom layer, the rock salt. 

In a similar situation for viscous fluids, the inversion of density leads to the 
so-called Ray leigh- Taylor instability due to buoyancy effect of gravity. In the 
numerical simulation by the use of SLA method, we shall present the results 
confirming the existence of similar instability for viscoelastic solid bodies. 

Consider a body consisting of two layers of elastic and viscoelastic solids 
as shown in Fig. [TJ The body is under the action of gravity g, and pi < p2- 

r, 



P2 


9 


Pi 





oooo 



oooo 



Figure 1: Boundary conditions and initial state. 
The upper boundary Fi is traction- free, the others T2 are roller-supported. The 
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initial state is an unstable equilibrium state, the movement of the salt-sediment 
interface can be initiated by a small perturbation. 

For illustrative purpose, we shall present numerical simulations in a two- 
dimensional domain for a Mooney-Rivlin type material. The proposed method 
has been applied to three-dimensional domain and some different class of vis- 
coelastic solid bodies with similar results. 

6.1 Formation of a salt diapir 

In this example we consider formation of a salt diapir initiated by a small 
perturbation at a very small region centered at the salt-sediment interface. 

• The dimension of the initial state is: length = 1,200 m, height of salt 
layer = 100 m, height of sediment layer = 200 m. 

• The material parameters for rock salt are: 

po = 2.2xl0 3 Kg/m 3 , si = 0, s 2 = -0.2xl0 3 Pa, 

A = -lO.OxlO 3 PaMa, fn = 15.0xl0 3 PaMa, /z 2 = M3 = 0. P = !0 9 Pa. 

• The material parameters for overburden sediment are: 

po = 3.0xl0 3 Kg/m 3 , si = 2.5xl0 3 Pa, s 2 = -7.5xl0 3 Pa, 
A = (J,i = M2 = M3 = 0, /3 = 10 9 Pa. 

• The incremental time: At = 0.1 Ma. 

The material data are of convenient choice for demonstration only. No attempt 
has been made to match the data to real properties of relevant materials. 

In Fig. [21 various stages of migration of the salt-sediment interface are shown, 
where n is the time step. One can easily see the formation of salt diapir as time 
step increases. The effect is primarily due to buoyancy force of density inversion. 
The diapir reaches its maximum height at about 20 Ma (million years) at n = 200 
with almost no change afterward as the diapir becomes mature. The formation 
of diapir is a result of very large deformation of the initial mesh as can be seen 
from Fig.[3]at n = 300. The deformed mesh is quite similar to the experimental 
results of silicone putty model of a diapir formed by spinning the model in a 
centrifuge by Dixon [TB]. 

6.2 Salt migration due to inclination 

In the second example, we consider two-layer structure of a greater extension, 
so that there are enough rock salt in the bottom layer to develop multi-diapirs. 
The migration is initiated by gradually lifting the base rock, which supports the 
salt layer, up to an inclination at an angle of one degree during the initial one 
million years (n = 10). 

• The dimension of the initial state is: length = 5,000 m, height of salt 
layer = 100 m, height of sediment layer = 200 m. 
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Figure 2: Migration of salt diapirism initiated by a small perturbation at the 
center of salt-sediment interface. The sequence represents the formation of a 
diapir during 100 million years (t — n At). 
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Figure 3: Original mesh and the deformed mesh at n — 300. 

• The material parameters for rock salt are: 

p = 2.2xl0 3 Kg/m 3 , Sl = 0, s 2 = -0.2xl0 3 Pa, 

A = -lO.OxlO 3 PaMa, fix = 15.0xl0 3 PaMa, \x 2 = fJ-3 = 0, /3 = 2xl0 9 Pa. 

• The material parameters for overburden sediment are: 

pa = 3.0xl0 3 Kg/m 3 , s x = 2.5xl0 3 Pa, s 2 = -7.5xl0 3 Pa, 
A = (J,i = Ma = M3 = 0, [3 = 2xl0 9 Pa. 

• The incremental time: At = 0.1 Ma. 

Due to the gravity, the lifting of the base rock on the left side pushes the 
body to the right which initiates the growth of a salt pillow at about the first 
50 million years (n = 500) as can be seen from Fig. |U As time goes on, an 
adjacent salt pillow appears as the first one becomes a salt diapir. The figures 
show the growth sequence from right to left of the appearance of different salt 
structures up to 150 million years (n — 1500). 

Such formation due to inclination of base rock was reported in permian salt 
complex of northern Germany as shown in the sketch (Fig.[S|) by Trusheim [17] . 
Although our numerical simulation is only two-dimensional, the similarity of the 
formation of salt structure is rather striking. However, since the numerical data 
are of convenient choice only, the estimated time in million of years in numerical 
simulation may not be of any real significance. 
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Figure 4: Salt migration initiated by a gradual lift of base rock to an angle of 
one degree at the initial 10 time steps on the left side. The figures show the 
growth sequence from right to the left for the appearance of salt structures. 
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Figure 5: Diagram of different types of salt structure due to a dip at an angle 
of more than one degree in Northern Germany (From F. Trushcim, Bull. Amer. 
Assoc. Petroleum Geologists 44 (I960)) 
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